8  Heteroscedasticidad

9 Introducción

Este documento presenta un análisis detallado sobre la heteroscedasticidad, uno de los problemas más comunes en el modelado de regresión que viola el supuesto de varianza constante de los errores. Se simulará un conjunto de datos que presenta heterocedasticidad, se mostrarán las técnicas gráficas y formales para su diagnóstico, y se explorarán cuatro soluciones principales para corregir o mitigar sus efectos: errores estándar robustos, Mínimos Cuadrados Ponderados (WLS), Mínimos Cuadrados Generalizados Factibles (FGLS) y la transformación de la variable de respuesta (Box-Cox).

9.0.1 Librerías Requeridas

Se cargan únicamente los paquetes necesarios para la ejecución de este análisis.

library(lmtest)   # Para bptest()
library(sandwich) # Para errores estándar robustos (vcovHC)
library(MASS)     # Para boxcox() y studres()

10 Ejemplo con Heteroscedasticidad

Se simula un conjunto de datos donde la varianza del término de error depende de una de las variables predictoras (\(Var(\epsilon_i | x_{i1}) = \sigma^2 x_{i1}\)), induciendo así heterocedasticidad.

num <- 1000
set.seed(109)
x1 <- rpois(num, 12)
x2 <- rgamma(num, 2)
# El término de error depende de x1, creando heterocedasticidad
error_term <- sapply(x1, function(i) rnorm(1, 0, sqrt(i))) 

# Ecuación del modelo
y <- 2*x1 + 0.5*x2 + error_term

# Se toma una muestra de los datos simulados
dat <- cbind(y, x1, x2)
dat <- as.data.frame(dat[sample(1:num, 500, replace = TRUE),])

11 Estimación Ordinaria de MCO

Se ajusta un modelo de regresión lineal estándar (Mínimos Cuadrados Ordinarios, MCO), ignorando inicialmente el problema de heterocedasticidad.

plot(dat)

reg <- lm(y ~ x1 + x2, data = dat)
summary(reg)

Call:
lm(formula = y ~ x1 + x2, data = dat)

Residuals:
     Min       1Q   Median       3Q      Max 
-10.3362  -2.1906   0.0218   2.1618  12.6354 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept) -0.41185    0.62178  -0.662    0.508    
x1           2.06254    0.04605  44.790  < 2e-16 ***
x2           0.45062    0.11348   3.971 8.22e-05 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 3.716 on 497 degrees of freedom
Multiple R-squared:  0.8032,    Adjusted R-squared:  0.8024 
F-statistic:  1014 on 2 and 497 DF,  p-value: < 2.2e-16

12 Diagnóstico de Homoscedasticidad

Se realizan pruebas gráficas y formales para detectar la presencia de heterocedasticidad.

• \(H_0\): Hay homocedasticidad (varianza de los errores constante).
• \(H_1\): Hay heterocedasticidad (varianza no constante).

# Diagnósticos gráficos
par(mfrow = c(2, 2))
plot(reg, which = 3) # Scale-Location plot
plot(dat$x1, reg$residuals, xlab="x1", ylab="Residuales") # Patrón de abanico
plot(dat$x2, reg$residuals, xlab="x2", ylab="Residuales")
plot(reg$fitted.values, reg$residuals, xlab="Valores Ajustados", ylab="Residuales")

par(mfrow = c(1, 1))

# Prueba formal: Test de Breusch-Pagan
bptest(reg)

    studentized Breusch-Pagan test

data:  reg
BP = 24.29, df = 2, p-value = 5.314e-06

Conclusión: Los gráficos (especialmente el de residuales vs. x1) muestran un patrón de abanico, donde la dispersión de los residuales aumenta con el valor del predictor. El Test de Breusch-Pagan confirma esto con un p-valor muy pequeño, rechazando la hipótesis nula de homocedasticidad. Se concluye que hay un problema de heteroscedasticidad.

---

13 Soluciones a la Heteroscedasticidad

13.1 Alternativa 1: Errores Estándar Robustos

Esta solución no cambia los coeficientes estimados (\(\hat{\beta}_j\)), pero corrige sus errores estándar para que la inferencia (p-valores, intervalos de confianza) sea válida a pesar de la heterocedasticidad. Se utiliza la matriz de varianza-covarianza de White (HC).

# Se recalculan las pruebas de hipótesis de los coeficientes con errores robustos
coeftest(reg, vcov. = vcovHC, type = "HC3")

t test of coefficients:

             Estimate Std. Error t value  Pr(>|t|)    
(Intercept) -0.411847   0.591242 -0.6966    0.4864    
x1           2.062544   0.048701 42.3515 < 2.2e-16 ***
x2           0.450623   0.109196  4.1267 4.314e-05 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

# Intervalos de confianza robustos
coefci(reg, vcov. = vcovHC, type = "HC3")

                2.5 %    97.5 %
(Intercept) -1.573489 0.7497948
x1           1.966860 2.1582284
x2           0.236080 0.6651654

# Comparación con los intervalos de confianza originales (no válidos)
confint(reg)

                 2.5 %    97.5 %
(Intercept) -1.6334920 0.8097976
x1           1.9720696 2.1530184
x2           0.2276652 0.6735802

# Comparación de modelos (equivalente a anova) con errores robustos
reg2 <- lm(y ~ x1, data = dat)
waldtest(reg2, reg, vcov = vcovHC)

Wald test

Model 1: y ~ x1
Model 2: y ~ x1 + x2
  Res.Df Df     F    Pr(>F)    
1    498                       
2    497  1 17.03 4.314e-05 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

13.2 Alternativa 2: Mínimos Cuadrados Ponderados (WLS)

Si se conoce la forma de la heterocedasticidad, se puede corregir el modelo ajustando por Mínimos Cuadrados Ponderados (WLS). En este caso, como \(Var(\epsilon_i | x_{i1}) \propto x_{i1}\), los pesos a utilizar son \(w_i = 1/x_{i1}\). WLS produce estimadores eficientes (BLUE).

# Se ajusta el modelo WLS especificando los pesos
reg_wls <- lm(y ~ x1 + x2, data = dat, weights = 1/x1)
summary(reg_wls)

Call:
lm(formula = y ~ x1 + x2, data = dat, weights = 1/x1)

Weighted Residuals:
    Min      1Q  Median      3Q     Max 
-3.6381 -0.6822  0.0115  0.6439  3.1569 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept) -0.55468    0.52680  -1.053    0.293    
x1           2.06977    0.04277  48.391  < 2e-16 ***
x2           0.47725    0.10813   4.414 1.25e-05 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 1.059 on 497 degrees of freedom
Multiple R-squared:  0.8278,    Adjusted R-squared:  0.8271 
F-statistic:  1194 on 2 and 497 DF,  p-value: < 2.2e-16

# Se comprueba que los residuales ponderados ahora son homocedásticos
plot(dat$x1, residuals(reg_wls) / sqrt(dat$x1), main = "Residuales Ponderados vs. x1")

plot(fitted(reg_wls), residuals(reg_wls) / sqrt(dat$x1), main = "Residuales Ponderados vs. Ajustados")

13.3 Alternativa 3: Mínimos Cuadrados Generalizados Factibles (FGLS)

Cuando la forma de la heterocedasticidad no se conoce, se puede estimar a partir de los residuales del modelo MCO inicial.

# Modelo FGLS simple (no recomendado): usa los residuales^2 como estimación de la varianza
fgls <- lm(y ~ x1 + x2, data = dat, weights = 1/reg$residuals^2) 
summary(fgls)

Call:
lm(formula = y ~ x1 + x2, data = dat, weights = 1/reg$residuals^2)

Weighted Residuals:
    Min      1Q  Median      3Q     Max 
-1.0566 -0.9986  0.8527  1.0012  2.1070 

Coefficients:
             Estimate Std. Error t value Pr(>|t|)    
(Intercept) -0.413320   0.034227  -12.08   <2e-16 ***
x1           2.063895   0.003033  680.51   <2e-16 ***
x2           0.440630   0.008257   53.37   <2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 0.9982 on 497 degrees of freedom
Multiple R-squared:  0.999, Adjusted R-squared:  0.9989 
F-statistic: 2.366e+05 on 2 and 497 DF,  p-value: < 2.2e-16

# Recomendación de Wooldridge: modelar el logaritmo de los residuales al cuadrado
modvar <- lm(log(reg$residuals^2) ~ x1 + x2, data = dat)
# Los pesos se obtienen como la inversa del exponencial de los valores ajustados de este modelo
fgls2 <- lm(y ~ x1 + x2, data = dat, weights = 1/exp(modvar$fitted.values))
summary(fgls2)

Call:
lm(formula = y ~ x1 + x2, data = dat, weights = 1/exp(modvar$fitted.values))

Weighted Residuals:
    Min      1Q  Median      3Q     Max 
-7.8342 -1.3529  0.0473  1.3365  5.7495 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept) -0.70667    0.53457  -1.322    0.187    
x1           2.07638    0.04523  45.903  < 2e-16 ***
x2           0.51789    0.10963   4.724 3.02e-06 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 2.089 on 497 degrees of freedom
Multiple R-squared:  0.8123,    Adjusted R-squared:  0.8115 
F-statistic:  1075 on 2 and 497 DF,  p-value: < 2.2e-16

# Se comprueba la homocedasticidad de los nuevos residuales ponderados
plot(dat$x1, studres(fgls2), main="Residuales Estudentizados del Modelo FGLS (Wooldridge)")

13.4 Alternativa 4: Transformar la Variable Respuesta (Box-Cox)

Otra estrategia es transformar la variable \(Y\) para estabilizar la varianza. El método de Box-Cox ayuda a encontrar la transformación de potencia óptima (\(\lambda\)).

\[ Y^{(\lambda)} = \frac{Y^\lambda - 1}{\lambda} \]

# La función boxcox() busca el valor de lambda que maximiza la verosimilitud
BC <- boxcox(reg)

lam <- BC$x[which.max(BC$y)]
lam

[1] 0.7878788

# Se transforma Y y se ajusta un nuevo modelo
yl <- (dat$y^lam - 1) / lam
dat2 <- cbind(dat, yl)
reg3 <- lm(yl ~ x1 + x2, data = dat2)
summary(reg3)

Call:
lm(formula = yl ~ x1 + x2, data = dat2)

Residuals:
    Min      1Q  Median      3Q     Max 
-6.2098 -1.1732  0.0298  1.1323  5.8288 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)  1.73006    0.31342   5.520 5.47e-08 ***
x1           1.04326    0.02321  44.945  < 2e-16 ***
x2           0.24004    0.05720   4.196 3.21e-05 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 1.873 on 497 degrees of freedom
Multiple R-squared:  0.8045,    Adjusted R-squared:  0.8037 
F-statistic:  1022 on 2 and 497 DF,  p-value: < 2.2e-16

# Se verifica si el nuevo modelo es homocedástico
plot(reg3, which = 3)

bptest(reg3)

    studentized Breusch-Pagan test

data:  reg3
BP = 10.27, df = 2, p-value = 0.005886

# Se compara la normalidad de los residuales del modelo original y el transformado
par(mfrow=c(1,2))
hist(studres(reg), main="Residuales Originales")
hist(studres(reg3), main="Residuales del Modelo Transformado")

par(mfrow=c(1,1))
shapiro.test(studres(reg3))

    Shapiro-Wilk normality test

data:  studres(reg3)
W = 0.99171, p-value = 0.006831

Conclusión: La transformación Box-Cox ha corregido exitosamente el problema de heterocedasticidad, como lo confirma el p-valor alto del test de Breusch-Pagan en el nuevo modelo.